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1 Summary 

Type la supernovae, i.e. stellar explosions which do not have hydrogen in their spectra, 
but intermediate-mass elements such as silicon, calcium, cobalt, and iron, have recently re- 
ceived considerable attention because it appears that they can be used as "standard candles" 
to measure cosmic distances out to billions of light years away from us. Observations of 
type la supernovae seem to indicate that we are living in a universe that started to accelerate 
its expansion when it was about half its present age. These conclusions rest primarily on 
phenomenological models which, however, lack proper theoretical understanding, mainly 
because the explosion process, initiated by thermonuclear fusion of carbon and oxygen into 
heavier elements, is difficult to simulate even on supercomputers. 

Here, we investigate a new way of modeling turbulent thermonuclear deflagration fronts 
in white dwarfs undergoing a type la supernova explosion. Our approach is based on a level 
set method which treats the front as a mathematical discontinuity and allows for full cou- 
pling between the front geometry and the flow field. New results of the method applied to 
the problem of type la supernovae are obtained. It is shown that in 2-D with high spatial 
resolution and a physically motivated subgrid scale model for the nuclear flames numeri- 
cally "converged" results can be obtained, but for most initial conditions the stars do not 
explode. In contrast, simulations in 3-D, do give the desired explosions and many of their 
properties, such as the explosion energies, lightcurves and nucleosynthesis products, are in 
very good agreement with observed type la supernovae. 



2 Introduction 

Numerical simulations of any kind of turbulent combustion have always been a challenge, 
and thermonuclear supernova explosions are no exception to that rule. This is mainly be- 
cause of the large range of length scales involved. In type la supernovae (SNe la), in particu- 
lar, the length scales of relevant physical processes range from 10~ 3 cm for the Kolmogorov- 
scale to several 10 7 cm for typical convective motions. In the currently favored scenario the 
explosion starts as a deflagration near the center of the star. Rayleigh-Taylor unstable blobs 
of hot burnt material are thought to rise and to lead to shear-induced turbulence at their 
interface with the unburnt gas. This turbulence increases the effective surface area of the 
flamelets and, thereby, the rate of fuel consumption; the hope is that finally a fast deflagra- 
tion might result, in agreement with phenomenological models of type la explosions. 

Despite considerable progress in the field of modeling turbulent combustion for astro- 
physical flows the correct numerical representation of the thermonuclear deflagration front 
has always been a weakness of the simulations. Methods used until recently were based 
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on the reactive-diffusive flame model, which artificially stretches the burning region over 
several grid zones to ensure an isotropic flame propagation speed. However, the soft tran- 
sition from fuel to ashes stabilizes the front against hydrodynamical instabilities on small 
length scales, which in turn results in an underestimation of the flame surface area and - 
consequently - of the total energy generation rate. Moreover, because nuclear fusion rates 
depend on temperature nearly exponentially, one cannot use the zone-averaged values of 
the temperature obtained this way to calculate the reaction kinetics. 

The front tracking method used in this project cures most of these weaknesses. It is 
based on the so-called level set technique which was originally introduced by Osher & 
Sethian (1988). They used the zero level set of a n-dimensional scalar function to represent 
(n - l)-dimensional front geometries. Equations for the time evolution of such a level set 
which is passively advected by a flow field are given in Sussman et al. (1994). The method 
has been extended to allow the tracking of fronts propagating normal to themselves, e.g. 
deflagrations and detonations (Smiljanovski et al., 1997). In contrast to the artificial broad- 
ening of the flame in the reaction-diffusion-approach, this algorithm is able to treat the front 
as an exact hydrodynamical discontinuity. 



3 The level set method 

The central aspect of our front tracking method is the association of the front geometry (a 
time-dependent set of points T) with an isoline of a so-called level set function G: 

r := {r | G(r) = 0). (1) 

Since G is not completely determined by this equation, we can additionally postulate that G 
be negative in the unburnt and positive in the burnt regions, and that G be a "smooth" func- 
tion, which is convenient from a numerical point of view. This smoothness can be achieved, 
for example, by the additional constraint that | VG| = 1 in the whole computational domain, 
with the exception of possible extrema and kinks of G. The ensemble of these conditions 
produces a G which is a signed distance function, i.e. the absolute value of G at any point 
equals the minimal front distance. The normal vector to the front is defined to point towards 
the unburnt material. 

The task is now to find an equation for the temporal evolution of G such that the zero 
level set of G behaves exactly as the flame. Such an expression can be obtained by the 
consideration that the total velocity of the front consists of two independent contributions: 
it is advected by the fluid motions at a speed v and it propagates normal to itself with a 
burning speed s. 

Since for deflagration waves a velocity jump usually occurs between the pre-front and 
post-front states, we must explicitly specify which state v and s refer to; traditionally, the 
values for the unburnt state are chosen. Therefore, one obtains for the total front motion 

D/ = v„ + s u ii. (2) 

The total temporal derivative of G at a point P attached to the front must vanish, since G is, 
by definition, always at the front: 

dG P dG dG 

^f = ir + VG -^ = ^ +D /- VG = < 3 > 

dt dt dt 

This leads to the desired differential equation describing the time evolution of G: 

dG 

ft-=- D '- VG - (4) 

This equation, however, cannot be applied on the whole computational domain, mainly 
because using this equation everywhere will in most cases destroy G's distance function 
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property. Therefore additional measures must be taken in the regions away from the front 
to ensure a "well-behaved" |VG| (Reinecke et al., 1999). 

The situation is further complicated by the fact that the quantities v„ and s u which 
are needed to determine D/ are not readily available in the cells cut by the front. In a 
finite volume context, these cells contain a mixture of pre- and post-front states instead. 
Nevertheless one can assume that the conserved quantities (mass, momentum and total 
energy) of the mixed state satisfy the following conditions: 



Here a denotes the volume fraction of the cell occupied by the unburnt state. In order to 
reconstruct the states before and behind the flame, a nonlinear system consisting of the 
equations above, the Rankine-Hugoniot jump conditions and a burning rate law must be 
solved. 




Fig. 1. Illustration of the basic principles of the level set method according to Smiljanovski et al. 
(1997): The piecewise linear front cuts the mixed cells into burnt and unburnt parts, a is the unburnt 
volume fraction of a cell, f} is the unburnt area fraction of a cell interface. The fluxes F„ and F b are 
calculated from the reconstructed states. 

Having obtained the reconstructed pre- and post-front states in the mixed cells, it is 
not only possible to determine D/, but also to separately calculate the fluxes of burnt and 
unburnt material over the cell interfaces. Consequently, the total flux over an interface can 
be expressed as a linear combination of burnt and unburnt fluxes weighted by the unburnt 
interface area fraction j3 (see Fig. 1): 



4 Implementation 

For our calculations, the front tracking algorithm was implemented as an additional module 
for the hydrodynamics code PROMETHEUS (Fryxell et al., 1989). Here we describe a 
simple implementation of most of the ideas discussed in the previous section which we 
will call "passive implementation". It assumes that the G-function is advected by the fluid 
motions and by burning and is only used to determine the source terms for the reactive Euler 
equations. It must be noted that there exists no real discontinuity between fuel and ashes in 
this case; the transition is smeared out over about three grid cells by the hydro-dynamical 
scheme, and the level set only indicates where the thin flame front should be. However, the 
numerical flame is still considerably thinner than in the reaction-diffusion approach. 

A complete implementation contains in-cell-reconstruction and flux-splitting as pro- 
posed by Smiljanovski et al. (1997) and outlined above. Therefore it describes exactly the 
coupling between the flame and the hydrodynamic flow. This generalized version of the 
code has been applied to hydrogen combustion in air by Reinecke et al. (1999) and, very 
recently, also to thermonuclear fusion by Ropke et al. (2003). 



p = ap u +(1 
pv = ap u \ u + (1 
pe = ap u e u +(1 



a)pb 
a)p b e b 



(5) 
(6) 
(7) 



F= J SF„ + (1 -J3)F b . 



(8) 
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4.1 G-Transport 

Since the front motion consists of two distinct contributions, it is appropriate to use an 
operator splitting approach for the time evolution of G. The advection term due to the fluid 
velocity vf reads in conservation form 

X^ 3r+ £ _VFpGjf=o (9) 

(Mulder et al., 1992). This equation is identical to the advection equation of a passive 
scalar, like the concentration of an inert chemical species. Consequently, this contribution 
to the front propagation can be calculated by PROMETHEUS directly. The additional flame 
propagation due to burning is calculated at the end of each time step and a re-initialization 
of G is done in order to keep it a signed distance function (see Reinecke et al. 1999). 



4.2 Source terms 

After the update of the level set function in each time step, the change of chemical compo- 
sition and total energy due to burning is calculated in the cells cut by the front. In order to 
obtain these values, the volume fraction a occupied by the unburnt material is determined 
in those cells by the following approach: from the value G,y and the two steepest gradients 
of G towards the front in x- and (/-direction a first-order approximation G of the level set 
function is calculated; then the area fraction of cell ij where G < can be found easily. 
Based on these results, the new concentrations of fuel, ashes and energy are obtained: 

X 'Ashes = max( 1 - QT, X As hes) ( 1 0) 

^Fuel = 1 ~ ^Ashes (H) 
Koi = e tot + 4(X A shes _ X Ashes) ( 1 2) 

In principle this means that all fuel found behind the front is converted to ashes and the 
appropriate amount of energy is released. The maximum operator in eq. (10) ensures that 
no "reverse burning" (i.e. conversion from ashes to fuel) takes place in the cases where 
the average ash concentration is higher than the burnt volume fraction; such a situation 
can occur in a few rare cases because of unavoidable discretization errors of the numerical 
scheme. 



5 Turbulent nuclear burning 

The system of equations described so-far can be solved provided the normal velocity of 
the burning front is known everywhere and at all times. In our computations it is deter- 
mined according to a flame-brush model of Niemeyer & Hillebrandt (1995a), which we 
will briefly outline for convenience. 

As was mentioned before, nuclear burning in degenerate dense matter is believed to 
propagate on microscopic scales as a conductive flame, wrinkled and stretched by local 
turbulence, but with essentially the laminar velocity. Due to the very high Reynolds num- 
bers macroscopic flows are highly turbulent and they interact with the flame, in principle 
down to the Kolmogorov scale. This means that all kinds of hydrodynamic instabilities feed 
energy into a turbulent cascade, including the buoyancy-driven Rayleigh-Taylor instability 
and the shear-driven Kelvin-Helmholtz instability. Consequently, the picture that emerges is 
more that of a "flame brush" spread over the entire turbulent regime rather than a wrinkled 
flame surface. For such a flame brush, the relevant minimum length scale is the so-called 
Gibson scale, defined as the lower bound for the curvature radius of flame wrinkles caused 
by turbulent stress. Thus, if the thermal diffusion scale is much smaller than the Gibson 
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scale (which is the case for the physical conditions of interest here) small segments of the 
flame surface are unaffected by large scale turbulence and behave as unperturbed laminar 
flames ("flamelets"). On the other hand side, since the Gibson scale is, at high densities, 
several orders of magnitude smaller than the integral scale set by the Rayleigh-Taylor ed- 
dies and many orders of magnitude larger than the thermal diffusion scale, both transport 
and burning times are determined by the eddy turnover times, and the effective velocity of 
the burning front is independent of the laminar burning velocity. 

A numerical realization of this general concept is presented in Niemeyer & Hillebrandt 
(1995a). The basic assumption was that wherever one finds turbulence this turbulence is 
fully developed and homogeneous, i.e. the turbulent velocity fluctuations on a length scale / 
are given by the Kolmogorov law u(l) = v(L)(l/L)^ 3 , where L is the integral scale, assumed 
to be equal to the Rayleigh-Taylor scale. Following the ideas outlined above, one can also 
assume that the thickness of the turbulent flame brush on the scale I is of the order of 
/ itself. With these two assumptions and the definition of the Gibson scale one finds for 

'gibs ^RT 

v(l) - s t (l) - it (/ gibs )(-U (13) 

V 'gibs / 

and d t (l) - I, where u('gibs) = ^lam defines Z g jb s , si am is the laminar burning speed and s t (l) is 
the turbulent flame velocity on the scale I. 

In a second step this model of turbulent combustion is coupled to our finite volume 
hydro scheme. Since in every finite volume scheme scales smaller than the grid size cannot 
be resolved, we express Z gibs in terms of the grid size A, the (unresolved) turbulent kinetic 
energy per unit mass, k sgs , and the laminar burning velocity: 

I 



- a[ 1 I (14) 



Here fe sgs is determined from a subgrid scale model (Clement, 1993; Niemeyer & Hille- 
brandt, 1995a) and, finally, the effective turbulent velocity of the flame brush on scale A is 
given by 

s t (A) = max(si am , v(A), u RJ ), (15) 
with u(A) = sJ2k sgs and u RT oc -^gA, where g is the local gravitational acceleration. 



6 Application to the supernova problem 

Different series of simulations were performed to check the numerical reliability of the 
employed models and to compare two- and three-dimensional explosions. 

6.1 Resolution study 

A crucial test for the validity of the models for the unresolved scales (in this case the flame 
and subgrid models) is to check the dependence of integral quantities, like the total energy 
release of the explosion, on the numerical grid resolution. Ideally, there should be no such 
dependence, indicating that all effects on unresolved scales are accurately modeled. 

Figure 2 shows the energy evolution of a centrally ignited white dwarf. The only dif- 
ference between the simulations is the central grid resolution, which ranges from 2 • 10 6 cm 
(model c3_2d_128) down to 2.5 • 10 5 cm (model c3_2d_1024). Model c3_2d_128 is obvi- 
ously under-resolved, but the results of the other calculations are in good agreement, with 
exception of the last stages, where the flame enters strongly non-uniform regions of the 
grid. 

So far, this kind of parameter study could only be performed in two dimensions, be- 
cause of the prohibitive cost of very highly resolved 3D simulations. Nevertheless the re- 
sults suggest that a resolution of « 10 6 cm should yield acceptable accuracy also in three 
dimensions. 
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Fig. 2. Time evolution of the total energy for identical initial conditions at different resolutions. While 
model c3_2d_128 is clearly under-resolved, the other simulations agree very well, at least in the early 
and intermediate stages. 
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Fig. 3. Burning front geometry evolution of model c3_2d_256. 



6.2 Comparison of 2D and 3D simulations 

In order to investigate the fundamental differences between two- and three-dimensional 
simulations, a 2D and a 3D model with identical initial conditions and resolution was cal- 
culated. Figures 3 and 4 show snapshots of the flame geometry at various explosion stages; 
the energy evolution of both models is compared in figure 5. 

It is evident that both simulations evolve nearly identically during the first few tenths 
of a second, as was expected. This is a strong hint that no errors were introduced into the 
code during the enhancement of the numerical models to three dimensions. At later times, 
however, the 3D calculation develops instabilities in the azimuthal direction, which could 
not form in 2D because of the assumed axial symmetry. As a consequence the total burning 
surface and the energy generation rate is increased, resulting in a higher overall energy 
release. 
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Fig. 4. Burning front geometry evolution of model c3_3d_256. One ring on the axes corresponds to 
10 7 cm. 
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Fig. 5. Comparison of the energy evolution for identical initial conditions in two and three dimen- 
sions. 
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6.3 The effect of different initial conditions 

In our approach, the initial white dwarf model (composition, central density, and velocity 
structure), as well as assumptions about the location, size and shape of the flame surface 
as it first forms fully determine the simulation results. At present, we have not changed 
the properties of the white dwarf but we concentrate on variations of the latter. In this 
context, the simultaneous runaway at several different spots in the central region of the 
progenitor star is of particular interest, a plausible ignition scenario suggested by Garcia- 
Senz & Woosley (1995). 




Fig. 6. Left panel: Snapshots of the flame front for a scenario with 5 ignition spots in 3D. The fast 
merging between the leading and trailing bubbles and the rising of the entire burning region is clearly 
visible. One ring on the coordinate axes corresponds to 10 7 cm. Right panel: The same as before, but 
for a high resolution model with 9 ignition spots. 

One simulation was carried out on a grid of 256 3 cells with a central resolution of 
10 6 cm and contained five bubbles with a radius of 3-10 6 cm, which were distributed ran- 
domly in the simulated octant within 1.6T0 7 cm of the star's center. In an attempt to re- 
duce the initially burned mass as much as possible without sacrificing too much flame 
surface, a very highly resolved second model was constructed. It contains nine randomly 
distributed, non-overlapping bubbles with a radius of 2T0 6 cm within 1.6T0 7 cm of the 
white dwarf's center. To properly represent these very small bubbles, the cell size was re- 
duced to A =5- 10 5 cm, so that a total grid size of 5 12 3 cells was required. 

During the first 0.5 seconds, the three models are nearly indistinguishable as far as 
the total energy is concerned (see Fig. 7), which at first glance appears somewhat surpris- 
ing, given the quite different initial conditions. A closer look at the energy generation rate 
actually reveals noticeable differences in the intensity of thermonuclear burning for the 
simulations, but since the total flame surface is initially very small, these differences have 
no visible impact on the integrated curve in the early stages. 

However, after about 0.5 seconds, when fast energy generation sets in, the nine-bubble 
model burns more vigorously due to its larger surface and therefore reaches a higher final 
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Fig. 7. Energy evolution of these three-dimensional explosion models (dashed and dashed-dotted). 
For comparison we also show the centrally ignited ("three fingers") model of Reinecke et al. (2002) 
(solid line) discussed earlier. 

energy level. Fig. 7 also shows that the centrally ignited model (c3_3d_256) is almost 
identical to the off-center model b5_3d_256 with regard to the explosion energetics. But, 
obviously, the scatter in the final energies due to different initial conditions appears to 
be small. Moreover, all models explode with an explosion energy in the range of what is 
observed. 



7 Predictions for observable quantities 

tn this Section we present a few preliminary results for various quantities which could, in 
principle, be observed and which therefore can serve as tests for the models. 

7.1 Lightcurves 

The most direct test of explosion models is provided by observed lightcurves and spectra. 
According to "Arnett's Law" lightcurves measure mostly the amount and spatial distribu- 
tion of radioactive 56 Ni in type la supernovae, and spectra measure the chemical composi- 
tion in real and velocity space. 

Sorokina & Blinnikov (2003) have used the results of one of our centrally ignited 3D- 
model, averaged over spherical shells, to compute colour lightcurves in the UBVt-bands. 
Their code assumes LTE radiation transport and loses reliability at later times (about 4 
weeks after maximum) when the supernova enters the nebular phase. Also, this assumption 
and the fact that the opacity is not well determined at longer wavelength make I-lightcurves 
less accurate. Keeping this in mind, the lightcurves shown in Fig. 8 look very promising. 
The main reason for the good agreement between the model and SN 1994D is the presence 
of radioactive Ni in outer layers of the supernova model at high velocities which is not be 
predicted by spherical models. 

7.2 Elemental and isotopic abundances 

A summary of the abundances obtained for all 3D models is given in Table 1 . Here "Mg" 
(as in Fig. 8) stands for intermediate-mass nuclei, and "Ni" for the iron-group, tn addition, 
the total energy liberated by nuclear burning is given. Since the binding energy of the white 
dwarf was about 5-10 50 erg, all models do explode. Typically one expects that around 80% 
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Fig. 8. UB Vl-colour lightcurves predicted by a centrally ignited 3D model (solid lines) in comparison 
with obseved data for a bright (98bu), a "normal" (94D), and a subluminous (91bg) supernova. 

of iron-group nuclei are originally present as 56 Ni bringing our results well into the range of 
observed Ni-masses. This success of the models was obtained without introducing any non- 
physical parameters, but just on the basis of a physical and numerical model of subsonic 
turbulent combustion. We also stress that our models give clear evidence that the often 
postulated deflagration-detonation transition is not needed to produce sufficiently powerful 
explosions. 



model name 


m Mg [M B ] 


n?Ni [M ] 


£ nl , [10 50 erg] 


c3_3d_256 


0.177 


0.526 


9.76 


b5_3d_256 


0.180 


0.506 


9.47 


b9_3d_512 


0.190 


0.616 


11.26 



Table 1. Overview over element production and energy release of all discussed supernova simulations 

Finally, we have "post-processed" one of our models in order to see whether or not also 
reasonable isotopic abundances are obtained. The results, shown in Fig. 9, are preliminary 
and should be considered with care. However, it is obvious that, with a few exceptions, also 
isotopic abundances are within the expected range and do not differ too much from those 
computed by means of phenomenological models, i.e. "W7". Exceptions include the high 
abundance of (unburned) C and O, and the overproduction of 48,50 Ti, 54 Fe, and 58 Ni. 

8 Supplementary studies 

In order to validate our approach to model the effects on unresolved scales, two studies were 
carried out. One tests several subgrid scale (SGS) models under the physical conditions 
encountered in SN la explosions and the other concerns the flame stability on small scales, 
where the turbulent cascade does not dominate the flame evolution. 
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Fig. 9. Isotopic abundances obtained for off-center ignited 3D model. 



8.1 Large eddy simulations of turbulent deflagration 

Subgrid scale models were tested in large eddy simulations (LES) of thermonuclear burn- 
ing in a turbulent flow. In these simulations, turbulence is driven by an artificial stochastic 
force field and the computational domain is cubic with periodic boundary conditions. The 
term large eddy simulation indicates that not all dynamic scales are numerically resolved. 
Compared to simulations of supernova explosions, there are two additional differences: 
Firstly, the periodic boundary conditions enforce a constant average mass density. Hence, 
there is no explosion in this setting. The burning process proceeds in the form of a percola- 
tion process as gradually all fuel is consumed and the system approaches a state in which 
degenerate C+O matter is replaced by 56 Ni and a particles in statistical nuclear equilibrium. 
Secondly, gravity is not included. Since the integral length scale is of the order 10 s cm, this 
is a sensible approximation. In consequence, we focus on the interaction of the deflagration 
with purely hydrodynamical turbulence for which the similarity theory of Kolmogorov ap- 
plies. Whether SGS models which do not account for effects induced by gravity are applica- 
ble to thermonuclear supernova explosions is still controversial. However, there are simple 
scaling arguments in favor of this conjecture. The single most important result emerging 
from our studies is that the evolution of the burning process is significantly affected by the 
SGS model in use. This underlines the importance of choosing a faithful SGS model and 
certainly bears consequences on the simulation of deflagration in SNe la. 

Our approach of determining the turbulent flame speed as outlined in Section 5 is based 
upon a dynamical equation for the SGS turbulence energy k sgs = \ql gs which can be casted 
into the following equation for the turbulence velocity q sgs : 

D 1 / \ 2 „, (I CA <??gs 

5^sgs- -V-(p4tfsgsVtfsgsJ-4|Vtfsgsl =(y\S \ - I - + — I q sgs d - — (16) 

The operator ^ is the Lagrangian time derivative Jj + v V. In this equation, several heuristic 
approximations, so-called closures, are incorporated. Associated with these closures are 
the length scales ( v = C v A eS l y2, ( e = 2y/2A ef {/C c and l K = C K A tS j y2. C A accounts 
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Fig. 10. Evolution of turbulent thermonuclear burning in a cubic domain subject to periodic boundary 
conditions. The panels show 2D contour sections of the normalized specific total energy, e tot /c^. The 
mean mass density is p ~ 2.903- 10 s gcrrr 3 and the initial sound speed c ~ 6.595 TO 8 cm s~' . Turbu- 
lence is produced artificially by a stochastic solenoidal force field. The bright regions of high specific 
energy contain burned material. The flame front is indicated by the thin white lines corresponding to 
the zero level set. 

for pressure effects in a compressible fluid. The closure parameters C v , C e , C K and Ca 
are a priori unknown. In the case of stationary isotropic turbulence in an incompressible 
fluid, values can be derived from analytic theories of turbulence. Particular examples are 
C v ~ 0.054, C e * 1.0 and C K ~ 0.1 (cf. Sagaut (2001), Section 4.3) and C A ~ -0.2 (Fureby 
et al., 1997). A e s is the effective length scale of the finite-volume scheme, in this case, the 
piece-wise parabolic method. Usually, this scale is set equal to the size A of the grid cells. 
However, from the energy spectra of direct numerical simulations (DNS) of forced isotropic 
turbulence, we concluded that is more appropriate to set A^ = f3A, where /3 is in the range 
1 .6 ... 1 .8 depending on the Mach number. The factor /3 accounts for the smoothing effect 
of the numerical scheme on top of the discretization. 

We have also determined statistical values of the closure parameters from DNS data. For 
subsonic flows, C v - 0.06, C e - 0.48 and C K - 0.36 appear to be representative values. The 
outcome of a LES of turbulent deflagration in a cube of 216 3 cells with these parameters is 
illustrated in Figure 10. Burning is ignited in eight small spherical zones at the beginning of 
the simulation and the fluid is set into motion by stochastic forcing. Turbulence is produced 
on a time scale T = L/V, where L = 2.16 • 10 s cm is the integral length scale and V = 
l00si am « 9.78 • 10 7 cms _I is the characteristic velocity of the flow. One can think of L 
being the typical size of the largest vortices generated by the stochastic force field and T is 
the associated autocorrelation time or, figuratively, the turn-over time. Initially, the burning 
zones are expanding very slowly as the flame propagation speed mostly equals the laminar 
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Fig. 11. Comparison of the evolution of SGS turbulence and thermonuclear burning for LES with 
three different variants of the SGS turbulence energy model. In one simulation, Clement's law of the 
wall was used to calculate the closure parameters, the second LES was computed with the constant 
parameters and for the third simulation, a semi-localized model was applied. In the top panels, the 
ratio of the mean SGS turbulence velocity to the RMS velocity and the laminar burning speed, re- 
spectively, is shown as function of normalized time t/T. In the middle row of panels, the first three 
statistical moments of the mass-weighted SGS turbulence velocity and, in the bottom panels, the time 
evolution of fuel, burning products and the mean energy generation rate are plotted. 



speed. After roughly one turn-over time has elapsed, SGS turbulence becomes increasingly 
space filling and enhances the flame propagation speed. At the same time, the folding and 
stretching of the flame front by the resolved flow greatly increases the surface and, thus, 
the rate of fuel consumption. These two effects appreciably accelerate the burning process 
as one can see from the evolution of the energy shown in Figure 10. At time t > 2T, most 
of the material has already been burnt. 

Constant closure parameters are a sensible choice for fully developed isotropic turbu- 
lence. However, significant deviations are expected for transient, intermittent or anisotropic 
flows. The deflagration in the cube is certainly transient and, naturally, the physical condi- 
tions in the vicinity of the flame front are anisotropic. In the case of a supernova explosion, 
these qualities are even more pronounced. For this reason, a more sophisticated method of 
determining the closure parameters is called for. Clement (1993) found for stellar structure 
modeling ad hoc rules to adjust the parameters C v and C e in order to account for compress- 
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ibility in stratified media. Since steep density gradients affect turbulence more or less like 
a wall, we shall refer to the rules proposed by Clement as "law of the wall". This method 
has been applied to SGS modeling in the simulations of thermonuclear supernovae as well. 
In comparison to the SGS turbulence energy model with constant parameters, Clement's 
law of the wall produces markedly different results as one can seen from the LES statistics 
plotted in Figure 1 1 . The mean of q sgs increases initially very rapidly and then settles at an 
almost constant level. Correspondingly, the average rate of energy release per unit volume 
due to thermonuclear burning, (B), increases more gradually as opposed to steep rise and 
high peak found for the SGS model with constant parameters. 

An entirely different method utilizes the filtering approach introduced by Germano 
(1992). Smoothing the resolved velocity field with a filter of characteristic length larger 
than A e g and invoking self-similarity assumptions, equations can be derived which yield 
closure parameters varying in space and time corresponding to the local structure and evo- 
lution of the flow. (Piomelli, 1993; Liu et al., 1994; Ghosal et al., 1995; Meneveau & 
Katz, 2000). This so-called dynamical procedure works particularly well for the produc- 
tion parameter C v . In the case of dissipation and diffusion, however, the method appears to 
be inadequate because of deficiencies in the respective closures. Therefore, we adopted a 
semi-localized approach. Results from a LES using this method, are shown in the panels on 
the very right of Figure 1 1 . The differences compared to the LES with constant parameters 
are not striking but Clement's method is clearly dismissed. The prediction of the latter that 
SGS turbulence is produced right at the beginning in advance of small-scale structure de- 
veloping in the resolved flow is physically unreasonable. Apart from that, once turbulence 
has developed in a particular spatial region, it becomes increasingly isotropic on decreasing 
length scales even in more complex flows. Thus, the remarkable similarity of the evolution 
of the burning process in the LES with constant and dynamical parameters, respectively. 
However, the choice of appropriated values for the closure parameters is notoriously diffi- 
cult, particularly, when it comes to the supernova explosion scenario in which turbulence 
is generated due to buoyancy effects originating from the density contrast between burnt 
and unburnt material. In conclusion, the semi-localized model is likely to be the preferable 
SGS model for the application to SNe la simulations after all. 

8.2 Investigation of the cellular burning regime 

Below the Gibson scale, flame propagation is not dominated by the turbulent cascade, since 
here the flame burns faster through turbulent eddies than these can deform it. On those 
scales flame evolution is determined by the (hydrodynamical) Landau-Darrieus (LD) insta- 
bility (Landau, 1944) and its counteracting nonlinear stabilization (Zel'dovich, 1966). In 
case of terrestrial flames this stabilization leads to a cellular steady-state pattern of the flame 
front giving rise to the cellular burning regime. After Niemeyer & Hillebrandt (1995b) have 
shown by means of hydrodynamical simulation that the LD instability acts under conditions 
of SNe la, we could confirm that also the cellular stabilization of the flame front holds for 
thermonuclear flames in white dwarfs. This was achieved in a series of 2D simulations. In 
order to reproduce hydrodynamical effects like the LD instability it turned out to be essen- 
tial to apply the complete implementation described in Sect. 3 (see Ropke et al. 2003). The 
stability of this cellular flame shape was tested for various fuel densities and for interaction 
of the flame with vortical flows. Two examples are given in the following. 

Figure 12 shows the propagation of an initially perturbed flame into quiescent fuel. 
The initial perturbations grow and in the nonlinear regime the flame exhibits a cellular 
structure. For the given setup, however, the cellular pattern of the same wavelength as 
the initial perturbation is not a stable solution. The snapshots at later times (Fig. 12b,c) 
illustrate the "merging" of cells until the steady-state structure of a single domain-filling 
cusp has formed. In case of higher numerical resolution, this fundamental flame structure 
may be superposed by a short-wavelength cellular pattern, which is advected toward the 
cusp and does not lead to a break-up of the domain-filling cell. This evolution is well in 
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Fig. 12. Flame propagation into quiescent fuel at p„ = 5 x 10 7 gem -3 , resolution: 200 x 200 cells, 
snapshots taken at (a) 0, (b) 7.5, (c) 15, and (d) 30 growth times tld of a perturbation with A = 
3.2 X 10 4 cm. Color-coded is the vorticity of the flow field. (The figure is also included in the color 
plates of this volume.) 

agreement with theoretical predictions. Simulations with different fuel densities p„ led to 
similar results. In the range of p„ = 1 x 10 7 g cirT 3 ... 1 x 10 9 g cm" 3 no significant flame 
destabilization could be observed. 

Flame interaction with a vortical fuel field is shown in Fig. 13. This is motivated by relic 
turbulent motions from the eddy cascade around the Gibson scale and from pre-ignition 
convective motions in the white dwarf. A parameter study with different fuel densities 
and various strengths of the velocity fluctuations in the fuel flow led to the result that a 
vortical flow of sufficient strength can break up the cellular pattern of the flame. In this 
case, however, no drastic self-turbulization of the flame but rather a smooth adaptation of 
the flame structure to the imprinted flow was observed. 

These results corroborate the assumption of large scale models, that the generation of 
turbulence is dominated by large-scale effects. In the cellular burning regime, the flame 
surface is enlarged compared to a planar configuration. This leads to an additional accel- 
eration of the flame. As a result of our study we suggest an increased lower cut-off of the 
flame velocity (cf. eq. (15)) instead of Mi am depending on background turbulence. 
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Fig. 13. Flame propagation into vortical fuel atp„ = 5 x 10 7 gcirr 3 ; velocity fluctuations at the right 
boundary: u'/"iam = 2.5; resolution: 300 x 200 cells; snapshots taken at (a) 0s, (b) 8.0 x 10~ 3 s, (c) 
1.6 x 10~ 2 s, and (d) 4.8 x 10~ 2 s. Color-coded is the vorticity of the flow field. (The figure is also 
included in the color plates of this volume.) 

9 Conclusion 

In this project, we have provided a new method to model the physics of thermonuclear com- 
bustion in degenerate dense matter of white dwarfs stars consisting of carbon and oxygen. 
Because not all relevant length-scales of this problem can be numerically resolved a numer- 
ical model to describe deflagration fronts with a reaction zone much thinner then the cells 
of the computational grid was presented. This new approach was applied to the simulate 
thermonuclear supernova explosions of Mchan white dwarfs in 2 and 3 dimensions. 

An implicit assumption of our numerical model is that on the resolved scales the flows 
are turbulent allowing us to describe the physics on the unresolved scales by a subgrid scale 
model, in the spirit of large-eddy simulations. The results presented here indicate that for 
supernova simulations this assumption is satisfied if we compute 3D models and use at 
least a 256 3 grid. The reason is simply that more structure on small length scales in better 
resolved models increases the rate of fuel consumption locally, but the turbulent velocity 
fluctuations are then smaller on the grid scale, compensating for this gain. 

All models we have computed (differing only in the ignition conditions and the grid 
resolution) explode. The explosion energy and the Ni-masses are only moderately depen- 
dent on the way the nuclear flame is ignited making the explosions robust. However, since 
ignition is a stochastic process, the differences we find may even explain some of the spread 
in observed SN la's. 

Based on our models we can predict lightcurves, spectra, and abundances, and the first 
preliminary results look promising. The lightcurves seem to be in very good agreement with 
observations, and also the nuclear abundances of elements and their isotopes are found to 
be in the expected range. Of course, the next step is to compute a grid of models, with 
varying white dwarf properties, and to compare them with the increasing data base of well- 
observed type la supernovae. The hope is that this will give us a tool to understand their 
physics and, thus, get confidence in their use for cosmology. 
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